function [field] = setUpField_v9_0(key)
%% setUpField_v9_0
%  Version 9.0
%  Author: Adeyinka Lesi
%  Date: 3/22/21
%  Project: Tumor Growth, Logarithmic Continuum Form
%  setUpField is a helper function to distributionHelper that sets up
%  some variables and puts them in a struct. The aim is to clean up the code
%  key: struct; contains variables used to set up the field
%  field: struct; contains variables used to run simulation
%  [field] = setUpField_v9_0(key)

%% Version History
%  2.0: Implementing irregular mesh around initial tumor size so initial
%  tumor is always on a grid point
%  3.0: setting up to handle multiple initial tumors
%  3.1: finding the index closest to 2
%  4.0: implementing new boundary condition model where N1 is distinguished
%  from the other sizes - N1 is basically treated as discrete. Will also do
%  the work of enabling irregular mesh
%  4.1: changing how initial values are assigned
%  4.2: case number is now stored in key struct
%  5.0: implementing size cutoff - any tumor size below cutoff is now known
%  and the probability for it will be distributed in the sizes below the
%  cutoff. Got rid of grid matching option.
%  5.1: added indexOf function; changed way tumors are put on field
%  5.2: changed x values to be median values of the bin rather than lower
%  values
%  5.3: implemented death starter
%  6.0: want to adjust position of partition and center points
%  6.1: 1) partitions at integer points
%       2) function values on half-intervals
%       3) initial distribution as a Gaussian curve or 
%  7.0: using transform function directly
%  7.1: version that works with getKey_backtrack, removed detection limit 
%  from getClosestIndex
%  7.2: slight modification to use key.DIST_SELECTOR choose which values in
%  key.INITIAL_TUMOR_DIST are used
%  7.3: changed format for key.DIST_SELECTOR to allow initial tumors from
%  time 2 data
%  8.0: To be used with getTransformedDistribution_CTC_v1_0, which allows
%  explicit CTC calculation and reseeding (requiring three new parameters).
%  Adding INITIAL_CTC_CONC, init_ctc
%  9.0: adding more options on placing unknown distributions
%  3/22/22: noticed error in indexOf where HARD_TUMOR_SIZE_LIMIT was used
%  instead of NUMBER_OF_SIZE_INTERVALS as max index value

display(['> ' mfilename]);
field = struct();

% def: name for saving data
field.name = [key.SET_NAME '_case' num2str(key.CASE)];

% def: transform function
trans = key.TRANSFORM;

% def: maximum tumor size allowed by code
field.xmax = key.HARD_TUMOR_SIZE_LIMIT;
field.ymax = trans.x2y(field.xmax);

% def: minimum tumor size
field.xmin = 1;
field.xmin2 = 2;
field.ymin = trans.x2y(field.xmin);
field.ymin2 = trans.x2y(field.xmin2);

% def: spatial discretization
field.dy = (field.ymax-field.ymin2)/(key.NUMBER_OF_SIZE_INTERVALS-1);
    
% def: transformed spatial variable
field.y = [field.ymin field.ymin2:field.dy:field.ymax];

% def: untransformed spatial variable
field.x = trans.y2x(field.y);

% generate initial distubution
% def: init is the initial tumor size distribution
field.init = zeros(length(field.y),1);
if(isfield(key,'TIME_ZERO_DIST'))
    init_dist = key.TIME_ZERO_DIST;    
elseif(isfield(key,'DIST_SELECTOR'))
    init_dist = key.INITIAL_TUMOR_DIST(key.DIST_SELECTOR);
else
    init_dist = key.INITIAL_TUMOR_DIST;
end
init_x = key.TIME_ZERO_SIZES;
init_y = trans.x2y(init_x);
init_indices = zeros(1,length(init_dist));
% instanciate CTC conc
if(isfield(key,'INITIAL_CTC_CONC'))
    field.init_ctc = key.INITIAL_CTC_CONC;
else
    field.init_ctc = 0;
end

% the floor function is used it is best to round the numbers used for it to
% prevent errors. floor_sg is the decimal place to round off to
floor_sg = 10;
% incorporate detection limit
field.detection_limit = key.DETECTION_LIMIT;
lowy = trans.x2y(field.detection_limit);
field.detection_index = floor(round((lowy-field.ymin2)/field.dy,floor_sg))+2;

% def: dys (y interval lengths), since we are enabling irregular mesh
field.dys = field.y(2:end) - field.y(1:end-1);
field.dxs = field.x(2:end) - field.x(1:end-1);

% the initial data may not be known entirely accurately - will allow the
% data to be placed on grid as a distribution. Either a Gaussian or a
% uniform distribution
% def: place_dist has a zero value if a uniform distribution is used and a
% value of 1 if a Gaussian distribution is used
if(isfield(key,'PLACEMENT_DIST'))
    place_dist = key.PLACEMENT_DIST;
else
    place_dist = zeros(1,length(init_dist));
end
% def: place_var is the half-width/std of the distribution
if(isfield(key,'PLACEMENT_RADIUS'))
    place_rad = key.PLACEMENT_RADIUS;
else
    place_rad = zeros(1,length(init_dist));
end

% must find where sizes falls among y values
% 'grid matching' is no longer used

for i = 1:length(init_dist)
    if init_dist(i) ~= 0
        % s is negative if init_x is below detection limit
        s = getClosestIndex(init_x(i),init_y(i),field);
        init_indices(i) = s;
    end
end

% update field.init(s), taking detection limit into account
for i = 1:length(init_indices)
    s = init_indices(i);
    % want to set this up so that the initial number calculated by
    % numerical integration is the intial number of tumors desired
    % currently using midpoint rule
    
    if s > 0
        if(place_rad(i) == 0)
            % placed in only one partition
            field.init(s) = field.init(s) + init_dist(i)/field.dxs(s);
        else
            % the values is distributed across partitions depending on
            % std and type of distribution
            % x0_dis is the size at which there is a tumor in discrete case
            % this implies that int^{x0_dis+1}_{x0_dis}{p dx} = 1 so that
            % (1)p(x0_dis+0.5) ~ 1; this is why 0.5 is added
            x0_dis = init_x(i); % position of tumor
            x0 = x0_dis + 0.5; % value at which distribution is specified
            rad = place_rad(i);
            switch place_dist(i)
                case 0
                    % uniform distribution
                    ylow = trans.x2y(x0-rad);
                    yhigh = trans.x2y(x0+rad);
                    low = getClosestIndex(x0-rad,ylow,field);
                    low = max(1,low);
                    high = getClosestIndex(x0+rad,yhigh,field);
                    high = max(1,high);
                    cdf = cdf_uniform(x0,rad,field.x(low:high+1));
                    cdf = cdf/(cdf(end)-cdf(1));
                    for j = low:high
                        frac = cdf(j-low+2) - cdf(j-low+1);
                        field.init(j) = field.init(j) + init_dist(i)*frac/field.dxs(j);
                    end
                case 1
                    % Gaussian distribution
                    ylow = trans.x2y(x0-rad);
                    yhigh = trans.x2y(x0+rad);
                    low = getClosestIndex(x0-rad,ylow,field);
                    low = max(1,low);
                    high = getClosestIndex(x0+rad,yhigh,field);
                    high = max(1,high);
                    cdf = cdf_gauss(x0,rad/6,field.x(low:high+1));
                    % adjust for fact that only 6 sigma is used
                    cdf = cdf/(cdf(end)-cdf(1));
                    for j = low:high
                        frac = cdf(j-low+2) - cdf(j-low+1);
                        field.init(j) = field.init(j) + init_dist(i)*frac/field.dxs(j);
                    end
                case 2
                    % wedge distribution (not symetric)
                    ylow = trans.x2y(x0-rad);
                    yhigh = trans.x2y(x0+rad);
                    low = getClosestIndex(x0-rad,ylow,field);
                    low = max(1,low);
                    high = getClosestIndex(x0+rad,yhigh,field);
                    high = max(1,high);
                    cdf = cdf_triangle(x0,rad,field.x(low:high+1));
                    cdf = cdf/(cdf(end)-cdf(1));
                    for j = low:high
                        frac = cdf(j-low+2) - cdf(j-low+1);
                        field.init(j) = field.init(j) + init_dist(i)*frac/field.dxs(j);
                    end
                case 3
                    % poisson distribution
                    ylow = trans.x2y(x0-rad);
                    yhigh = trans.x2y(x0+rad);
                    low = getClosestIndex(x0-rad,ylow,field);
                    low = max(1,low);
                    high = getClosestIndex(x0+rad,yhigh,field);
                    high = max(1,high);
                    cdf = cdf_poisson(x0,field.x(low:high+1));
                    cdf = cdf/(cdf(end)-cdf(1));
                    for j = low:high
                        frac = cdf(j-low+2) - cdf(j-low+1);
                        field.init(j) = field.init(j) + init_dist(i)*frac/field.dxs(j);
                    end
            end
        end
        %     else
        %         % all indices below detection are assumed to have a given
        %         % distribution
        %         switch place_dist(i)
        %             case 0
        %                 % uniform distribution
        %                 low = 1;
        %                 high = field.detection_index-1;
        %                 rad = (field.x(high+1)-field.x(low))/2;
        %                 xm = (field.x(high+1)+field.x(low))/2;
        %                 cdf = cdf_uniform(xm,rad,field.x(low:high+1));
        %                 cdf = cdf/(cdf(end)-cdf(1));
        %                 for j = low:high
        %                     frac = cdf(j-low+2) - cdf(j-low+1);
        %                     field.init(j) = field.init(j) + init_dist(i)*frac/field.dxs(j);
        %                 end
        %             case 1
        %                 % gaussian distribution
        %                 low = 1;
        %                 high = field.detection_index-1;
        %                 rad = (field.x(high+1)-field.x(low))/12;
        %                 xm = (field.x(high+1)+field.x(low))/2;
        %                 cdf = cdf_gauss(xm,rad,field.x(low:high+1));
        %                 cdf = cdf/(cdf(end)-cdf(1));
        %                 for j = low:high
        %                     frac = cdf(j-low+2) - cdf(j-low+1);
        %                     field.init(j) = field.init(j) + init_dist(i)*frac/field.dxs(j);
        %                 end
        %             case 2
        %                 % declining right triangle distribution  (not symetric)
        %                 low = 1;
        %                 high = field.detection_index-1;
        %                 rad = (field.x(high+1)-field.x(low))/2;
        %                 xm = (field.x(high+1)+field.x(low))/2;
        %                 cdf = cdf_triangle(xm,rad,field.x(low:high+1));
        %                 cdf = cdf/(cdf(end)-cdf(1));
        %                 for j = low:high
        %                     frac = cdf(j-low+2) - cdf(j-low+1);
        %                     field.init(j) = field.init(j) + init_dist(i)*frac/field.dxs(j);
        %                 end
        %             case 3
        %                 % poisson distribution
        %                 low = 1;
        %                 high = field.detection_index-1;
        %                 xl = field.x(low);
        %                 xm = 18+xl+6*sqrt(9+xl);
        % %                 xr = field.x(high+1)-field.x(low);
        % %                 xm = 18+xr-6*sqrt(9+xr);
        %                 rad = sqrt(xm);
        %                 cdf = cdf_poisson(xm,rad,field.x(low:high+1));
        %                 cdf = cdf/(cdf(end)-cdf(1));
        %                 for j = low:high
        %                     frac = cdf(j-low+2) - cdf(j-low+1);
        %                     field.init(j) = field.init(j) + init_dist(i)*frac/field.dxs(j);
        %                 end
        %         end
    end
end

% def: init_indices contains the indices that have tumors initially
field.init_indices = consolidateWithNoZeros(init_indices);

% def: soft_limit is the index of the max size initially considered
highest_index = field.init_indices(end);
field.soft_limit = min(highest_index + key.SOFT_LIMIT_BUFFER, ...
                       length(field.y) - 2);
field.buffer = key.SOFT_LIMIT_BUFFER;
field.reset_threshold = key.LIMIT_RESET_THRESHOLD;

% calculate derivatives
% def: derivatives of transform
% field.dydx = [field.dys./field.dxs 0];
field.dydx = trans.dydx(field.x);
field.d2ydx2 = trans.d2ydx2(field.x);

% generate time parameters
% def: t_max is the time duration over which the tumor distribution will be
%   observed
field.tmax = key.TIME_LIMIT;
% def: dt is the time step used for the numerical simultion
field.dt = key.TIME_STEP;

% define surgery parameters
if(key.SURGERY_TIME < key.TIME_LIMIT)
    field.surg_t = key.SURGERY_TIME; % time surgery occurs
else
    field.surg_t = key.TIME_LIMIT;
end
field.surg_limit_x = key.SURGERY_CUTOFF_SIZE; % min size removable by surgery
field.surg_limit_y = trans.x2y(field.surg_limit_x);
field.surg_index = floor(round((field.surg_limit_y-field.ymin2)/field.dy,floor_sg))+2;
field.using_surgery = key.CONDUCTING_SURGERY;

% set up death_starter
field.using_death_starter = key.USING_DEATH_STARTER;
field.death_start_time = key.DEATH_START_TIME;

% setting up indexOf function for use externally
field.indexOf = @indexOf;

disp('-> End of setup');
end

function index = getClosestIndex(x_size,y_size,field)
% getClosestIndex gives the index of the y coordinate closest to the size
% must find where init_size falls among y values

if(x_size >= field.x(2))
    index = floor(round((y_size-field.y(2))/field.dy,10))+2;
else
    index = 1;
end

index = min(index,length(field.x)-1);
end

function index = indexOf(x,key)
% indexOf gives the index of the partition the value falls into
% x: vector, size quantity
% key: struct, contains all required data
y = key.TRANSFORM.x2y(x);
% want to round of sizes to nearest index
pn = (y-key.FIELD.y(2))/key.FIELD.dy;
index = floor(round(pn,10))+2;
index = min(index,key.NUMBER_OF_SIZE_INTERVALS);
index = max(index,1);
end

function pithy = consolidateWithNoZeros(redun)
vec = zeros(size(redun));
count = 0;
for i = 1:length(redun)
    if(redun(i) ~= 0 && ~ismember(redun(i),vec))
        count = count + 1;
        vec(count) = redun(i);
    end
end
pithy = vec(1:count);
end

% cdf for uniform distribution of width 2*r centered at x0
% % the distribution is symetric around x=1
function vals = cdf_uniform(x0,r,xlist)
    vals = zeros(1,length(xlist));
    for i = 1:length(xlist)
        if(xlist(i) < x0-r)
            vals(i) = 0;
        elseif(xlist(i) < x0+r)
            vals(i) = (xlist(i)-x0+r)/2/r;
        else
            vals(i) = 1;
        end
        if(xlist(i) < 2-x0-r)
            % no change
        elseif(xlist(i) < 2-x0+r)
            vals(i) = vals(i)+(xlist(i)+x0+r-2)/2/r;
        else
            vals(i) = vals(i)+1;
        end
    end
end

% cdf for Gaussian distribution of std r centered at x0
% the distribution is symetric around x=1
function vals = cdf_gauss(x0,r,xlist)
    vals = 0.5*(erf((xlist-x0)/r/sqrt(2))+erf((xlist+x0-2)/r/sqrt(2)));
end

% cdf for declining right angle triangular distribution of width  2*r
% centered at x0
function vals = cdf_triangle(x0,r,xlist)
    vals = (xlist-x0+r)/4/r^2.*(3*r-xlist+x0);
end

% cdf for poisson distribution centered at x0
function vals = cdf_poisson(x0,xlist)
    % round the numbers since floor function is used
    xlist = round(xlist*1e10)*1e-10;
    vals = poisscdf(xlist-xlist(1),x0-xlist(1));
end
